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Abstract. Instabilities and strong dynamical interactions between mul- 
tiple giant planets have been proposed as a possible explanation for the 
surprising orbital properties of extrasolar planetary systems. In particu- 
lar, dynamical instabilities seem to provide a natural mechanism for pro- 
ducing the highly eccentric orbits seen in many systems. Previously, we 
performed numerical integrations for the dynamical evolution of planetary 
systems containing two giant planets of equal masses initially in nearly 
circular orbits very close to the dynamical stability limit. We found the 
ratio of collisions to ejections in these simulations was greater than the 
ratio of circular orbits to eccentric orbits among the known extrasolar 
planets. Further, the mean eccentricity of the planets remaining after an 
ejection was larger than the mean eccentricity of the known extrasolar 
planets. Recently, we have performed additional integrations, general- 
izing to consider two planets of unequal masses. Our new simulations 
reveal that the two-planet scattering model can produce a distribution 
of eccentricities consistent with the observed eccentricity distribution for 
plausible mass distributions. Additionally, this model predicts a maxi- 
mum eccentricity of about 0.8, in agreement with observations. Early 
results from simulations of three equal-mass planets also reveal a reduced 
frequency of collisions and a broad range of final eccentricities for the 
retained inner planet. 



1. Background 

The known extrasolar planets (Fig. 1 left) can be roughly divided into two 
groups: those with short-period, nearly circular orbits (a <^ 0.07 AU) and those 
with wider and more eccentric orbits (a ^ 0.07 AU). Many of the short-period 
planets, like their prototype 51 Peg, are so close to their parent star that tidal 
dissipation would have likely circularized their orbits, even if they were originally 
eccentric (Rasio et al. 1996). Thus, their small observed eccentricities do not 
provide a good indicator of their dynamical history. However, the large eccen- 
tricities of most planets with longer periods require an explanation. A planet 
that would have formed from a protoplanetary disk in the standard manner is 
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unlikely to have developed such a large eccentricity, since dissipation in the disk 
tends to circularize orbits. Dynamical instabilities leading to the ejection of 
one planet while retaining another planet of comparable mass could explain the 
observed distribution of eccentricities. 

In previous papers (Rasio & Ford 1996; Ford, Havlickova, & Rasio 2001), we 
conducted Monte Carlo integrations of planetary systems containing two equal 
mass planets initially placed just inside the Hill stability limit (Gladman 1993). 
We numerically integrated the orbits of such systems until there was a collision, 
or one planet was ejected from the system, or some maximum integration time 
was reached. The two most common outcomes were collisions between the two 
planets, producing a more massive planet in a nearly circular orbit between the 
two initial orbits, and ejections of one planet from the system while the other 
planet remains in a tighter orbit with a large eccentricity. The relative frequency 
of these two outcomes depends on the ratio of the planet radius to the initial 
semi-major axis. While this model could naturally explain how planets might 
acquire large eccentricities, upon comparing our results to the observed planets 
we found two important differences. First, for the relevant radii and semi- 
major axes, we found that collisions are more frequent in our simulations than 
nearly circular orbits among the presently known extrasolar planets. Second, 
of the systems which ended with one planet being ejected from the system, the 
eccentricity distribution of the remaining planet was concentrated in a narrow 
range which is greater than the typical eccentricity of the known extrasolar 
planets (See Fig. 1, right). Here we present results for new simulations involving 
two planets with unequal masses and three planets with equal masses. 



2. Two Planets, Unequal Masses 
2.1. Numerical Setup 

The new simulations use the mixed variable symplectic algorithm (Wisdom & 
Holman 1991) modified to allow for close encounters between planets as imple- 
mented in the publicly available code Mercury (Chambers 1999). The results 
presented below are based on ~ 10 4 numerical integrations. 

Throughout the integrations, close encounters between any two bodies were 
logged, allowing us to present results for any values of the planetary radii using 
a single set of orbital integrations. We consider a range of radii to allow for the 
uncertainty in both the physical radius and the effective collision radius allowing 
for dissipation in the planets. When two planets collided, mass and momentum 
conservation were used to compute the final orbit of the resulting single planet. 

Each run was terminated when one of the following four conditions was 
encountered: (i) one of the two planets became unbound (which we defined as 
having a radial distance from the star of 1000 a^init); (ii) a collision between 
the two planets occurred assuming Ri/a\^ n \t = ^ m in/oq,init = 0.1i?j up /5AU = 
0.95 x 10~ 5 ; (iii) a close encounter occurred between a planet and the star 
(defined by having a planet come within r m \ n / ai,\mt = IOJSq/IAU = 0.06 of 
the star); (iv) the integration time reached t max = 5 • 10 6 — 2 • 10 7 depending on 
the masses of the planets. These four types will be referred to as "collisions," 
meaning a collision between the two planets, "ejections," meaning that one 
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Figure 1. Left: Orbital eccentricity versus semi- major axis of known 
extrasolar planets. The eccentricity of planets in very short period 
orbits is limited by tidal damping. Nearly all of the remaining plan- 
ets have large eccentricities. The area of the circle is proportional to 
msini. Planets represented by shaded circles are part of a multiple 
planet system. Right: Cumulative distribution of eccentricity. The 
solid line shows the eccentricity distribution of the remaining planet 
after the other planet has been ejected from the system in our previ- 
ous simulations with two planets of equal mass. The dotted line shows 
the eccentricity distribution of the known extrasolar planets (excluding 
those with P < lOd). Among systems resulting in an ejection, the final 
eccentricity is larger than typical for the observed extrasolar planets. 



planet was ejected to infinity, "star grazers," meaning that one planet had a 
close pericenter passage, and "two planets," meaning that two bound planets 
remained in a (possibly new) dynamically stable configuration. 

Our numerical integrations were performed for a system containing two 
planets, with mass ratios < nii/M < 10~ 2 , where rrii is the mass of one of 
the planets and M is the mass of the central star. A mass ratio of m/M ~ 10~ 3 
corresponds to m ~ lMj up for M = 1M@. The initial semimajor axis of 
the outer planet (a2,init) was drawn from a uniform distribution ranging from 
0.9 • a^init (1 + A c ) to a^init (1 + A c ), where 1 + A c is the critical ratio above 
which Hill stability is guaranteed for initially circular coplanar orbits (Gladman 
1993). The initial eccentricities were distributed uniformly in the range from 
to 0.05, and the initial relative inclination in the range from 0° to 2°. All 
remaining angles (longitudes and phases) were randomly chosen between and 
2-7T. Throughout this paper we quote numerical results in units such that G = 
a>i init = M = 1. In these units, the initial orbital period of the inner planet is 
Pi' ~ 2vr. 

2.2. Results 

Branching Ratios For the range of rrii/M considered (up to 10~ 2 ), the fre- 
quency of the various outcomes is relatively insensitive to the ratio of the total 
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planet mass to the stellar mass. However the frequencies of the different out- 
comes are significantly effected by the planetary mass ratio and by the planetary 
radii. Collisions are most common for large planetary radii and mass ratios near 
unity, while ejections are most common for small planetary radii and mass ratios 
far from unity. The frequency of outcomes does not depend on whether the more 
massive planet is initially in the inner or outer orbit. In Fig. 2 (left), we show 
the frequency of ejection or collision as a function of the planetary radii with 
solid and dashed lines, respectively. The black lines show the result for equal 
mass planets. The other colors are for different mass ratios. Ejections become 
more common and collisions less frequent as the planetary mass ratio departs 
from unity. 

Collisions Collisions leave a single, larger planet in orbit around the star. The 
energy in the center-of-mass frame of the two planets is much smaller than 
the binding energy of a giant planet. Therefore, we model the collisions as 
completely inelastic and assume that the two giant planets simply merge together 
while conserving total momentum and mass. Under this assumption, we have 
calculated the distributions of orbital parameters for the collision products. The 
final orbit has a semi-major axis between the two initial semi-major axes, a 
small eccentricity, and a small inclination. While collisions between planets may 
affect the masses of extrasolar planets, a single collision between two massive 
planets does not cause significant orbital migration or eccentricity growth if 
the planets are initially on low-eccentricity, low-inclination orbits near the Hill 
stability limit. 

Ejections When one planet is ejected from the system, the less massive planet 
is nearly always ejected if < 0.8, where m < and m > refer to the masses of the 
less and more massive planets, respectively. Since the escaping planet typically 
leaves the system with a very small (positive) energy, energy conservation sets 
the final semimajor axis of the remaining planet slightly less than 

£/ „ Q 2 m < (1) 

ai ai (mi + m 2 ) + (a 2 - a±) mi ' 

Thus, the ejection of one of two equal mass planets results in the most significant 
reduction in the semi-major axis, but is limited to a ^ f > 0.5. The remaining 

planet acquires a significant eccentricity, but its inclination typically remains 
small. The eccentricity distribution for the remaining planet is not sensitive 
to the sum of the planet masses, but depends significantly on the mass ratio. 
In Fig. 2 (right) we show the cumulative distributions for the eccentricity after 
a collision for different mass ratios in different colors. While any one mass 
ratio results in a narrow range of eccentricities, a distribution of mass ratios 
would result in a broader distribution of final eccentricities. Also note that 
there is a maximum eccentricity which occurs for equal mass planets. Thus, 
the two planet scattering model predicts a maximum eccentricity of about 0.8 
independent of the distribution of planet masses. While this compares favorably 
with the presently known planets, future observations will certainly test this 
prediction. 
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Figure 2. Left: Branching ratio versus planet radius. The frequency 
of collisions (dashed lines) and ejections (solid lines) depends on the 
planet's radius relative to the semi-major axis (x-axis) and the plane- 
tary mass ratio (different colors) . The fraction of integrations resulting 
in one planet being ejected and the other planet remaining in an ec- 
centric orbit is increased when the two planetary masses differ. Right: 
Cumulative distribution of eccentricity. Each line shows the cumula- 
tive distribution of the eccentricity of the remaining planet after the 
other planet has been ejected from the system in our simulations for a 
particular mass ratio. The distributions are not sensitive to the total 
mass of the planets or to which planet is more massive. 

Stargrazers In a small fraction of our numerical integrations (~ 3%) one planet 
underwent a close encounter with the central star. Due to the limitations of the 
numerical integrator used, the accuracy of our integrations for the subsequent 
evolution of these systems can not be guaranteed. Moreover, these systems 
could be affected by additional forces (e.g., tidal forces, interaction with the 
quadrupole moment of the star) that are not included in our simulations and 
would depend on the initial separation and the radius of the star. Nevertheless, 
our simulations suggest that for giant planets with initial semimajor axes of a few 
AU, the extremely close pericenter distances necessary for tidal circularization 
around a main sequence star (leading to the formation of a 51-Peg-type system) 
are possible, but rare. 

2.3. Discussion 

We now consider whether the the two planet scattering model could produce 
a distribution of eccentricities consistent with that of known extrasolar plan- 
ets with periods longer than 10 days. In Fig. 3 (left) we show the cumula- 
tive distribution for the eccentricity of the remaining planet in our simulations. 
Here we have chosen a simple distribution of planet masses, P(m) ~ m _1 for 
lAfjup < m < 10Mj U p, which is consistent with the mass distribution of the 
observed extrasolar planets (Tabachnik & Tremaine 2001). The red dashed line 
is just for systems which resulted in an ejection, while the green dotted line 



6 



Ford, Rasio, <fe Yu 



includes systems which resulted in a collision. Both are reasonably close to the 
observed distribution shown by the solid black line. It should be noted that the 
mass distribution used for this model has two cutoffs. While the 10Mj up cutoff 
is supported by observations, the lower cutoff is not constrained by present ob- 
servations and is unlikely to be physical. Alternative mass distributions differing 
for m < lMj U p could alter this result. A detailed comparison to observations 
would require the careful consideration of observational selection effects, as well 
as the unknown initial mass distributions. Nevertheless, we find that the two- 
planet scattering model is able to reproduce the eccentricity distribution of the 
known extrasolar planets for plausible mass distributions. 




Figure 3. Left: Simulation with mass distribution. Here we com- 
pare the observed distribution of eccentricities (excluding planets with 
P < 10 d, solid line) to the eccentricity distributions produced by 
our simulations when each planet's mass is drawn independently from 
P(m) ~ m _1 for lMj up < m < 10Mj up . The dashed line shows 
the eccentricity of the planet remaining after an ejection, while the 
dotted line includes systems resulting in a collision. Right: Results 
from simulations for 3 equal-mass planets. Here we compare the ob- 
served distribution of eccentricities (excluding those with P < 10 d, 
solid black line) to the eccentricity distribution of the remaining inner 
planet after there has been one ejection or collision for one of two sets 
of simulations with different initial eccentricities. The dashed blue and 
dotted red lines are for simulations with e^t = 0. and e; n it = 0.05, 
respectively. 



3. Three Planets, Equal Masses 

We have begun performing simulations of unstable systems containing three 
equal-mass planets using the conservative Bulirsch-Stoer integrator in Mercury 
(Chambers 1999). The results presented here are based on ~ 10 3 numerical 
integrations. For three-planet systems, we set the radius of each planet at the 
beginning of the integration. For determining when collisions occurred, each 
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planet was assigned a radius of Ri = i?j up . A planet was removed from the sys- 
tem if its separation from the central star became greater than 1000AU (ejected) 
or less than r m - m = R & (collided with star). 

For three-planet systems, a single ejection or collision results in a system 
containing two planets. The remaining two planets may either be dynamically 
stable or undergo a subsequent ejection or collision to end with a single planet. 
Each run was terminated when two of the planets were removed by being ejected, 
colliding with another planet, or colliding with the star, or when the maximum 
integration time was reached t max = 2 x 10 6 years. 

Our numerical integrations were performed for a system containing three 
planets, each with a mass ratio m/M ~ 10 -3 . The initial semi-major axes are 
fixed at 1, 1.39, and 1.92 AU, so that the ratio of orbital periods is P3/P2 = 
P2/P1 = 1-63. Each planet is initially placed on a nearly circular orbit with an 
eccentricity of either 0.0 or 0.05 and an inclination of ±1°. All remaining angles 
(longitudes and phases) were randomly chosen between and 2tt. 

The typical time for instability to develop is ~ 10 6 years. In the majority 
of cases, two planets were left in a stable configuration after one planet was 
ejected or two planets collided. In Fig. 3 (right) we show the distributions of the 
eccentricity of the inner planet after a collision or ejection of one planet. It is in 
reasonable agreement with the observed distribution for extrasolar planets. Both 
the eccentricity distribution (Fig. 3, right) and the corresponding distribution 
of the inclination of the inner planet (Fig. 4, left) are in qualitative agreement 
with those presented by Marzari & Weidenschilling (2002, based on slightly 
different initial conditions). In Fig. 4 (right) we show the cumulative distribution 
for Aw = |(f&2 + 0J2) — ('tli + the angle between the periapses of the two 
planets remaining after one planet has been ejected from the system. We do 
not find the two remaining planets to have preferentially aligned periapses (cf. 
Malhotra 2002). 

4. Conclusions 

Dynamical interactions between planets of unequal mass reduces the frequency of 
collisions as compared to scattering between equal-mass planets. In particular, 
planet-planet scattering can reproduce the observed distribution of eccentric 
orbits using a plausible distribution of mass ratios. To determine if collisions 
are sufficiently rare will require a more careful comparison to the observations. 
Independent of the distribution of planet masses, planet-planet scattering by 
two planets initially on nearly circular orbits produces a maximum eccentricity 
of ~ 0.8. While this compares favorably with the presently known extrasolar 
planets, future detections will test this prediction of the two-planet scattering 
model. Alternatively, the observed properties of extrasolar planets could also 
be produced as a result of dynamical instabilities developing in three-planet 
systems. 
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Figure 4. Left: Cumulative distribution of inclination. Each line 
shows the distribution of the inclination of the remaining inner planet 
(with respect to the initial orbital plane) after there has been one ejec- 
tion or collision for one of two sets of simulations with different initial 
eccentricities. The dashed blue and dotted red lines are for simula- 
tions with ei n it = 0. and e- m \t = 0.05, respectively.. Right: Cumulative 
distribution of the periapse misalignment angle between two retained 
planets. 
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